subroutine writeinfo(fname, app, info)
use molecule
use iounits
implicit none
   !----------------------------------------------------------------------------
   character (len=*)   :: fname, app
   integer, intent(in) :: info
   !----------------------------------------------------------------------------
   integer             :: i, j, k, l, idum, idum2, id, ioerr
   character(len=100)  :: fmtstr
   !----------------------------------------------------------------------------
   !
   if (info.ne.6) then
      fname = trim(fname)//".dat"
   else
      fname = trim(fname)//".xyz"
   endif

   open(ioout,file=trim(fname),position=trim(app),action='write',iostat=ioerr)
   !
   select case (info)
   case ( 1 )
      call writeatom()
   case ( 2 )
      call writebond()
   case ( 3 )
      call writeangle()
   case ( 4 )
      call writedihed()
   case ( 5 )
      call WriteAll()
   case ( 6 )
      call WriteXyz()
   end select
   !
   close(ioout)
   return
contains
   !
   subroutine writeatom() ! Atomic info
   implicit none
      !-------------------------------------------------------------------------
      if (natom.lt.1.or.natomtype.lt.1) return
      idum2 = 0
      do i = 1, natom
         if ( typeid(types(i)).gt.0 ) idum2 = idum2 + 1
      enddo
      !
      idum = int(log10(dble(idum2)+0.1))+1
      fmtstr = '(/,10x,I??,5x,"atoms")'
      write(fmtstr(9:10), '(I2.2)') idum
      write(ioout, fmtstr) idum2
      !
      idum = int(log10(dble(unique(typeid))+0.1))+1
      fmtstr = '(/,10x,I??,5x,"atom types",/)'
      write(fmtstr(9:10), '(I2.2)') idum
      write(ioout, fmtstr) unique(typeid)
      fmtstr = '("#",2x,I4,2x,A2)'
      do i = 1, natomtype
         if ( typeid(i).gt.0 ) write(ioout, fmtstr) typeid(i), typename(i)
      enddo
      write(ioout,'(/,"# atom type in sequence: ",/,"# ",$)')
      fmtstr = '(I??,", ",$)'
      idum = int(log10(dble(natomtype)))+1
      write(fmtstr(3:4),'(I2)') idum
      do i = 1, natom
         write(ioout, fmtstr) typeid( types(i) )
      enddo
      !
      if ( allocated( molpos ) ) then
         fmtstr = '(2(I??,1x),3(1x,F20.10), "  # ", A2)'
         idum   = int(log10(dble(natom)))+1
         write(fmtstr(5:6),'(I2.2)') idum
         write(ioout, '(/,"# LAMMPS atomic style position")')
         idum = 0
         do i = 1, natom
            id = typeid( types(i) )
            if ( id.gt.0 ) then
               idum = idum + 1
               write(ioout, fmtstr) idum, id, molpos(:,i), typename(types(i))
            endif
         enddo
         !
         write(ioout, '(/,"# To code into alkene")')
         fmtstr = '("atposMol(:,",I??,")= (/",F10.5,"D0,",F10.5,"D0,",F10.5,"D0 /)")'
         write(fmtstr(17:18),'(I2.2)') int(log10(dble(natom)))+1
         idum = 0
         do i = 1, natom
            id = typeid( types(i) )
            if ( id.gt.0 ) then
               idum = idum + 1
               write(ioout, fmtstr) idum, molpos(:,i)
            endif
         enddo
         !
         write(ioout, '(/)')
      endif
   return
   end subroutine writeatom

   subroutine writebond()
   implicit none
      !-------------------------------------------------------------------------
      if (nbondtype.lt.1.or.nbond.lt.1) return
      idum = int(log10(dble(unique(bondid))+0.1))+1
      fmtstr = '(/,10x,I??,5x,"bond types",/)'
      write(fmtstr(9:10), '(I2.2)') idum
      write(ioout, fmtstr) unique(bondid)
      fmtstr = '("#",2x,I4,4x,A2,"-",A2)'
      do i = 1, nbondtype
         if ( bondid(i).gt.0 ) write(ioout, fmtstr) bondid(i), bondname(1:2,i)
      enddo
      write(ioout,'(//,"Bonds",/)')
      !
      fmtstr = '(4I??)'
      idum = int(log10(dble(max(nbond,natom))))+2
      write(fmtstr(4:5),'(I2.2)') idum
      idum = 0
      do i = 1, nbond
         id = bondid(bonds(3,i))
         if ( id.gt.0 ) then
            idum = idum + 1
            write(ioout, fmtstr) idum, id, bonds(1:2,i)
         endif
      enddo
      idum2 = int(log10(dble(idum)))+1
      fmtstr = '(/,10x,I??,5x,"bonds",/)'
      write(fmtstr(9:10),'(I2.2)') idum2
      write(ioout,fmtstr) idum
      !
      write(ioout,'(/,"# Bond list to code into alkene")')
      fmtstr = '("bondsMol(:,",I??,") = (/",I??,2(", ",I??), "/)" )'
      idum = int(log10(dble(max(nbond,natom))))+2
      write(fmtstr(17:18),'(I2.2)') idum
      write(fmtstr(30:31),'(I2.2)') idum
      write(fmtstr(41:42),'(I2.2)') idum
      idum = 0
      do i = 1, nbond
         id = bondid(bonds(3,i))
         if ( id.gt.0 ) then
            idum = idum + 1
            write(ioout, fmtstr) idum, id, bonds(1:2,i)
         endif
      enddo
      !
   return
   end subroutine writebond

   subroutine writeangle()
   implicit none
      !-------------------------------------------------------------------------
      if (nangletype.lt.1.or.nangle.lt.1) return
      idum = int(log10(dble(unique(angleid))+0.1))+1
      fmtstr = '(/,10x,I??,5x,"angle types",/)'
      write(fmtstr(9:10), '(I2.2)') idum
      write(ioout, fmtstr) unique(angleid)
      fmtstr = '("#",2x,I4,4x,A2,"-",A2,"-",A2)'
      do i = 1, nangletype
         if ( angleid(i).gt.0 ) write(ioout, fmtstr) angleid(i), anglename(1:3,i)
      enddo
      write(ioout,'(//,"Angles",/)')
      !
      fmtstr = '(5I??)'
      idum = int(log10(dble(max(nangle,natom))))+2
      write(fmtstr(4:5),'(I2.2)') idum
      idum = 0
      do i = 1, nangle
         id = angleid(angles(4,i))
         if ( id.gt.0 ) then
            idum = idum + 1
            write(ioout, fmtstr) idum, id, angles(1:3,i)
         endif
      enddo
      idum2 = int(log10(dble(idum)))+1
      fmtstr = '(/,10x,I??,5x,"angles",/)'
      write(fmtstr(9:10),'(I2.2)') idum2
      write(ioout,fmtstr) idum
      !
      write(ioout,'(//,"# Bond list to code into alkene")')
      fmtstr = '("AnglsMol(:,",I??,") = (/",I??,3(", ",I??)," /)")'
      idum = int(log10(dble(max(nangle,natom))))+2
      write(fmtstr(17:18),'(I2.2)') idum
      write(fmtstr(30:31),'(I2.2)') idum
      write(fmtstr(41:42),'(I2.2)') idum
           idum = 0
      do i = 1, nangle
         id = angleid(angles(4,i))
         if ( id.gt.0 ) then
            idum = idum + 1
            write(ioout, fmtstr) idum, id, angles(1:3,i)
         endif
      enddo
      !
   return
   end subroutine writeangle

   subroutine writedihed()
   implicit none
      !-------------------------------------------------------------------------
      if (ndihedtype.lt.1.or.ndihed.lt.1) return
      idum = int(log10(dble(unique(dihedid))+0.1))+1
      fmtstr = '(/,10x,I??,5x,"dihedral types",/)'
      write(fmtstr(9:10), '(I2.2)') idum
      write(ioout, fmtstr) unique(dihedid)
      fmtstr = '("#",2x,I4,4x,A2,"-",A2,"-",A2,"-",A2)'
      do i = 1, ndihedtype
         if ( dihedid(i).gt.0 ) write(ioout, fmtstr) dihedid(i), dihedname(1:4,i)
      enddo
      write(ioout,'(//,"Dihedrals",/)')
      !
      fmtstr = '(6I??)'
      idum = int(log10(dble(max(ndihed,natom))))+2
      write(fmtstr(4:5),'(I2.2)') idum
      idum = 0
      do i = 1, ndihed
         id = dihedid(diheds(5,i))
         if ( id.gt.0 ) then
            idum = idum + 1
            write(ioout, fmtstr) idum, id, diheds(1:4,i)
         endif
      enddo
      idum2 = int(log10(dble(idum)))+1
      fmtstr = '(/,10x,I??,5x,"dihedrals",/)'
      write(fmtstr(9:10),'(I2.2)') idum2
      write(ioout,fmtstr) idum
      !
      fmtstr = '("DihrsMol(:,",I??," ) = (/",I??,4(", ",I??)," /)")'
      idum = int(log10(dble(max(ndihed,natom))))+2
      write(fmtstr(17:18),'(I2.2)') idum
      write(fmtstr(31:32),'(I2.2)') idum
      write(fmtstr(42:43),'(I2.2)') idum
      idum = 0
      do i = 1, ndihed
         id = dihedid(diheds(5,i))
         if ( id.gt.0 ) then
            idum = idum + 1
            write(ioout, fmtstr) idum, id, diheds(1:4,i)
         endif
      enddo
   return
   end subroutine writedihed

   ! To get the number of unique non-zero elements in a 1D integer array
   integer function unique( ids )
   implicit none
      !-------------------------------------------------------------------------
      integer, intent(in)  :: ids(:)
      integer, allocatable :: id_uniq(:)
      integer              :: num, it, jt, id
      logical              :: found
      !-------------------------------------------------------------------------
      unique = 0
      num = size(ids)
      if ( allocated( id_uniq ) ) deallocate( id_uniq )
      allocate( id_uniq(num) )
      do jt = 1, num
         id = ids(jt)
         if (id.eq.0) cycle
         found = .false.
         do it = 1, unique
            if (id_uniq(it).eq.id) found = .true.
         enddo
         if (.not.found) then
            unique = unique + 1
            id_uniq(unique) = id
         endif
      enddo
   end function unique

   subroutine WriteAll() ! Write all info that suits code alkene
   implicit none
      !-------------------------------------------------------------------------
      integer :: ntypeMol, nbndTypMol, nangtypMol, ndihtypMol
      integer :: natomMol, nbondsMol, nanglsMol, ndihrsMol
      !
      if (natom.lt.1.or.natomtype.lt.1) return
      write(ioout, '(A)') trim(molname)
      write(ioout, '("Molecular Formula")')

      ntypeMol   = unique(typeid)
      if (nbondtype.lt.1.or.nbond.lt.1) then
         nbndTypMol = 0
      else
         nbndTypMol = unique(bondid)
      endif
      if (nangletype.lt.1.or.nangle.lt.1) then
         nangtypMol = 0
      else
         nangtypMol = unique(angleid)
      endif
      if (ndihedtype.lt.1.or.ndihed.lt.1) then
         ndihtypMol = 0
      else
         ndihtypMol = unique(dihedid)
      endif
      idum = max(ntypeMol, nbndTypMol, nangtypMol, ndihtypMol)
      fmtstr = '(4(1x,I??))'
      write(fmtstr(8:9), '(I2.2)') int(log10(dble(idum)+0.1))+1
      write(ioout, fmtstr) ntypeMol, nbndTypMol, nangtypMol, ndihtypMol
      natomMol = 0
      do i = 1, natom
         if ( typeid(types(i)).gt.0 ) natomMol = natomMol + 1
      enddo
      nbondsMol = 0
      do i = 1, nbond
         if ( bondid(bonds(3,i)).gt.0 ) nbondsMol = nbondsMol + 1
      enddo
      nanglsMol = 0
      do i = 1, nangle
         if ( angleid(angles(4,i)).gt.0) nanglsMol = nanglsMol + 1
      enddo
      ndihrsMol = 0
      do i = 1, ndihed
         if ( dihedid(diheds(5,i)).gt.0 ) ndihrsMol = ndihrsMol + 1
      enddo
      idum = max(natomMol, nbondsMol, nanglsMol, ndihrsMol)
      fmtstr = '(4(1x,I??))'
      write(fmtstr(8:9), '(I2.2)') int(log10(dble(idum)+0.1))+1
      write(ioout, fmtstr) natomMol, nbondsMol, nanglsMol, ndihrsMol
      !
      write(ioout, '("1.0  4")' ) ! default density and first dimension of torsion coeffs
      !
      write(ioout, '(/,"Elemental Info")')
      idum = 0
      fmtstr = '(I??,1x,A," mass sigma epsilon dist-min")'
      write(fmtstr(3:4),'(I2.2)') int(log10(dble(maxval(typeid))+0.1))+1
      do i = 1, natomtype
         if ( typeid(i).gt.0 ) then
            idum = idum + 1
            write(ioout,fmtstr) typeid(i), trim(typename(i))
         endif
      enddo

      if (.not.allocated(chg)) then
         allocate(chg(natom))
         chg = 0.D0
      endif
      write(ioout, '(/,"Atoms")')
      if ( allocated( molpos ) ) then
         fmtstr = '(3(I??,1x), F7.4, 3(1x,F12.8), "  # ", A2)'
         idum = max(maxval(typeid), natom)
         write(fmtstr(5:6),'(I2.2)') int(log10(dble(idum)))+1
         idum = 0
         do i = 1, natom
            id = typeid( types(i) )
            if ( id.gt.0 ) then
               idum = idum + 1
               write(ioout, fmtstr) idum, 1, id, chg(i), molpos(:,i), typename(types(i))
            endif
         enddo
      endif

      write(ioout, '(/,"Bond Coeffs")')
      fmtstr = '(I4,2x,"k_b  r_0  #", A2,"-",A2)'
      do i = 1, nbondtype
         if ( bondid(i).gt.0 ) write(ioout, fmtstr) bondid(i), bondname(1:2,i)
      enddo

      write(ioout, '(/,"Bonds")')
      fmtstr = '(4I??)'
      idum = int(log10(dble(max(nbond,natom))))+2
      write(fmtstr(4:5),'(I2.2)') idum
      idum = 0
      do i = 1, nbond
         id = bondid(bonds(3,i))
         if ( id.gt.0 ) then
            idum = idum + 1
            write(ioout, fmtstr) idum, id, bonds(1:2,i)
         endif
      enddo
      !
      write(ioout, '(/,"Angle Coeffs")')
      fmtstr = '(I4,2x,"k_theta theta_0  #", A2,"-",A2,"-",A2)'
      do i = 1, nangletype
         if ( angleid(i).gt.0 ) write(ioout, fmtstr) angleid(i), anglename(1:3,i)
      enddo
      write(ioout,'(//,"Angles")')
      fmtstr = '(5I??)'
      idum = int(log10(dble(max(nangle,natom))))+2
      write(fmtstr(4:5),'(I2.2)') idum
      idum = 0
      do i = 1, nangle
         id = angleid(angles(4,i))
         if ( id.gt.0 ) then
            idum = idum + 1
            write(ioout, fmtstr) idum, id, angles(1:3,i)
         endif
      enddo
      !
      write(ioout, '(/,"Dihedral Coeffs")')
      fmtstr = '(I4,2x,"Dihedral-coefficients  #", A2,"-",A2,"-",A2,"-",A2)'
      do i = 1, ndihedtype
         if ( dihedid(i).gt.0 ) write(ioout, fmtstr) dihedid(i), dihedname(1:4,i)
      enddo
      write(ioout,'(//,"Dihedrals")')
      fmtstr = '(6I??)'
      idum = int(log10(dble(max(ndihed,natom))))+2
      write(fmtstr(4:5),'(I2.2)') idum
      idum = 0
      do i = 1, ndihed
         id = dihedid(diheds(5,i))
         if ( id.gt.0 ) then
            idum = idum + 1
            write(ioout, fmtstr) idum, id, diheds(1:4,i)
         endif
      enddo
   return
   end subroutine WriteAll

   subroutine WriteXyz() ! Atomic info
   implicit none
      !-------------------------------------------------------------------------
      if (natom.lt.1.or.natomtype.lt.1) return
      idum2 = 0
      do i = 1, natom
         if ( typeid(types(i)).gt.0 ) idum2 = idum2 + 1
      enddo
      !
      write(ioout, '(I0)') idum2
      write(ioout, '(A)') trim(molname)
      !
      if ( allocated( molpos ) ) then
         fmtstr = '(A2,3(1x,F20.10))'
         idum = 0
         do i = 1, natom
            id = typeid( types(i) )
            if ( id.gt.0 ) then
               idum = idum + 1
               write(ioout, fmtstr) typename(types(i)),  molpos(:,i)
            endif
         enddo
      endif
   return
   end subroutine WriteXyz
end subroutine
